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Abstract 

We are interested in the fast computation of the exact value of integrals of polynomial 
functions over convex polyhedra. We present speed ups and extensions of the algorithms 
presented in previous work by some of the authors. We present a new software implemen- 
tation and provide benchmark computations. The computation of integrals of polynomials 
over polyhedral regions has many applications; here we demonstrate our algorithmic tools 
solving a challenge from combinatorial voting theory. 
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1. Introduction 

Integration is a fundamental problem arising in many contexts, from engineering to 
algebraic geometry. For example, center of mass, moments and products of inertia about 
various axes are integrals used in any type of dynamic simulation or modeling such as 
computer games [1, 2, 3]), similarly, normalized volumes indicate the degrees of toric vari- 
eties and closely related to moment maps of symplectic manifolds [4, 5]). Integration over 
polyhedra is particularly useful because many domains can be approximated by polyhedra 
and then decomposed into convex polyhedra (e.g., a tetrahedral mesh decomposition etc.). 

In this work we are interested in the exact evaluation of integrals over convex polyhedral 
regions. More precisely, let P be a <i-dimensional rational convex polyhedron in IR n and let 
/ G Q[xi, . . . ,x n ] be a polynomial with rational coefficients. We consider the problem of 
efficiently computing the exact value of the integral of the polynomial / over P, which we 
denote by j p f dm. Here we use the integral Lebesgue measure dm on the affine hull (P) 
of the polytope P. This general setting is quite useful because, when the polytope is full- 
dimensional, the integral Lebesgue measure coincides with the standard Riemann integral 
but generalizes it naturally to cases when the polytope is not full-dimensional. Another 
reason to compute in this setting is that the volume of P and every integral of a polynomial 
function with rational coefficients yields rational numbers. Finally this normalization of 
the measure occurs naturally in Euler-Maclaurin formulas for a polytope P, which relate 
sums over the lattice points of P with integrals over the various faces of P. We remark 
that the computer algebra community has dedicated a great deal of effort to develop a 
different kind of exact integration, understood to be the automatic computation of the 
antiderivatives of functions, as predicted by the fundamental theorem of Calculus [6], but 
we do not discuss this kind of exact integration here. 

Regarding the theoretical computational complexity of our problem, it is very educa- 
tional to look first at the case when / is the constant polynomial 1, and the answer is 
simply a volume. It has been proved that already computing the volume of polytopes of 
varying dimension is #P-hard [7, 8, 9, 10, 11], and that even approximating the volume 
is hard [12]. More recently in [13] it was proved that computing the center of mass of a 
polytope is #P-hard. 

We report on a new C++ implementation, LattE integrale [14], which extends the 
work done in [15] and [16]. This paper is mostly an experimental and practical study, 
but it also slightly develops the theory of [15]. This article presents useful formulas for 
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integration of power of linear forms over simplicial cones that complement those presented 
in [15]. 

Our method of computation relies on powerful mathematical ideas. It was proved in 
[15] that, unlike general polynomials, integrals over simplices of arbitrary powers of linear 
forms, or of polynomials with a fixed number of variables, can be computed in polynomial 
time. In this case our algorithms use known properties of integrals of exponentials of 
linear forms (see [17], [18]). This allows for fast calculation over general polytopes using 
two methods that depend on two different decompositions of polyhedra. General polyhedra 
can be decomposed as either a disjoint union of simplices, i.e., triangulations, or as signed 
cone decompositions of the kind proposed by Brion, Lasserre, Lawrence, and Varchenko 
[19, 20, 10, 21]. The polynomial-time complexity for integration over simplices shown in 
[15] can be extended to more polyhedra as long as their decompositions are of "small" size 
(note that this is always the case in fixed dimension). 

This paper is organized as follows. We begin in Section 2 recalling the mathematical 
ideas at the heart of our algorithms (although we omit details of proofs, they can be found in 
the references). We begin with a short review of polyhedral geometry, specially valuations. 
In Section 3 we discuss details about the implementation including main subroutines and 
data structures. In Section 4 we first discuss speed improvements for integrating over 
simplices from earlier work in [15] , and then we report on several benchmarks of integration 
over arbitrary polytopes. More experimental tables are available online [14]. We conclude 
our paper with an application: we solve a computational challenge from combinatorial 
voting theory. 

2. Mathematical preliminaries 

In this section we recall the necessary mathematical background used in our algorithms. 
We state results without proof but excellent background sources for what is going to be 
discussed here include [15, 22, 23, 24] and the references mentioned there. 

2.1. Polyhedra and polynomials 

A convex rational polyhedron P in M. d (we will simply say polyhedron) is defined as 
the intersection of a finite number m of closed half spaces bounded by rational affine 
hyperplanes. We say that P is full- dimensional (in M d ) if the affine span of P is Mr. 
When P = { x : Ax < b } for a m x d matrix A and m- vector b, P is said to be given 
by a halfspace or h-representation. When P is the convex hull conv(\^) of finitely many 
points in IR d , V = {vi, . . .v n }, P is said to be given by a vertex or v-representation. We 
can switch between the h- and v-representations of a <i-dimensional polyhedron using well- 
known algorithms (see [25, 26]). A polytope P is a compact polyhedron. A cone C is a 
polyhedral cone (with vertex 0) and an affine cone is a translation s + C of a cone C. 
A cone C is called simplicial if it is generated by linearly independent vectors of lR d . A 
simplicial cone C is called unimodular if it is generated by linearly independent integral 
vectors vi,...,Vk such that {vx, . . . , Vk} can be completed to an integral basis of Z d . An 
affine cone C is called simplicial (respectively, simplicial unimodular) if the associated cone 



3 



is. The set of vertices of P is denoted by V(P). For each vertex s G V(P), the cone of 
feasible directions C S (P) at the vertex s is the cone of all vector y such that s + ey G P for 
some e > 0. The tangent cone of a polytope P at a vertex s is the affine cone s + C S (P) 
(this is a translation of C S (P)). For details in all these notions see, e.g., [22]. 

For the integration of a full-dimensional polytope we consider the standard Lebesgue 
measure, which gives volume 1 to the fundamental domain of the lattice Z n . But if P lies 
inside an affine subspace L + a, with L a rational linear subspace of dimension n < d, we 
will normalize the Lebesgue measure on L, so that the volume of the fundamental domain 
of the lattice L PI Z d is 1. Thus for any affine subspace L + a parallel to L, we define the 
integral Lebesgue measure dm by translation. For example, the diagonal of the unit square 
has length 1 instead of \/2. This has the great advantage that for rational input our output 
will always be an (exact) rational number J p f dm in the usual binary encoding. 

One important point of our method is that all computations are done in the representa- 
tion polynomials given by powers of linear forms. It is well-known that any homogeneous 
polynomial of degree M can be decomposed as a sum of M-th powers of linear forms. For 
example, one can decompose the polynomial / as a sum / = ce(£, x) M with at most 2 M 
terms. This decomposition is given by the following well-known identity for monomials: If 
x* /L = x? 1 x?*---x% n , then 

* m = -L £ ( _ 1) im,-( P1+ ...^)^a.. w blXl + ... + PnXn) iM, (1) 

1 0< Pi <Mi VPn ^ 

where |M| = Mi + • • ■ + M n < M. Of course, when dealing with general polynomials, 
this same formula can be applied for as many monomials as is necessary. For example, the 
polynomial 7x 2 + y 2 + 5z 2 + 2xy + 9yz can be written as |(12(2x) 2 - 9(2y) 2 + (2z) 2 + 8(x + 
y) 2 + m{y + z) 2 ). 

It is worth noting that the above formula does not yield an optimal decomposition, but 
it suffices to generate a polynomial-time algorithm on fixed degree |M| or fixed number of 
variables [15]. The problem of finding a decomposition with the smallest possible number 
of summands is known as the polynomial Waring problem. What is the smallest integer 
r(M, n) such that a generic homogeneous polynomial f(x\, . . . , x n ) of degree Minn vari- 
ables is expressible as the sum of r(M, n) M-th powers of linear forms? This problem was 
solved for general polynomials by Alexander and Hirschowitz [27] (see [28] for an extensive 
survey), but there is no computational or constructive version of this result that would 
yield the optimal decomposition for an specific input polynomial and the bounds may be 
much too pessimistic on the average situation. Only very recently Carlini et al. [29] gave 
efficient decompositions of a monomial. However, their decomposition involves roots of 
unity, and here we are interested in an arithmetic version of the problem where everything 
is expressed using rational forms and rational coefficients. But we can see that the explicit 
rational construction we use in our code is not too far away from the optimum. 

Table 1 lists the average number of powers of linear forms necessary to decompose 
monomials of given degree generated uniformly at random. To create the monomials, we 
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Table 1: Average number of powers of linear forms plus or minus one standard deviation 
necessary to express one monomial in d variables, averaged over 50 monomials of the same 
degree 



Monomial Degree 



d 



10 



20 



30 



40 



50 



3 14 ±3 (6.6 ± 1.2) 

4 16 ±5 (1.1 ±0.2) 

5 19 ±4 (1.5 ±0.4) 

6 20 ±5 (2.0 ±0.6) 

7 21 ±5 (2.4 ±0.9) 

8 21 ±5 (2.9 ±0.9) 
10 24 ±5 (3.5 ±1.1) 



x 10 1 
x 10 2 
x 10 2 
x 10 2 
x 10 2 
x 10 2 
x 10 2 



4.0 ± 0.5) x 10 2 (1.2 ± 0.1) x 10 3 

1.1 ± 0.2) x 10 3 (4.5 ± 0.6) x 10 3 

2.2 ± 0.6) x 10 3 (1.2 ± 0.3) x 10 4 
4.1 ± 1.2) x 10 3 (3.2 ± 0.8) x 10 4 
6.7 ± 2.4) x 10 3 (7.1 ± 2.1) x 10 4 
1.1 ± 0.5) x 10 4 (1.4 ± 0.5) x 10 5 (9.8 ± 2.7) x 10 5 
2.1 ± 0.9) x 10 4 (4.1 ± 1.6) x 10 5 (4.5 ± 1.7) x 10 6 



(2.7 ±0.2) x 10 3 (5.2 ±0.2) x 10 3 
(1.3 ±0.2) x 10 4 (3.0 ±0.2) x 10 4 



(4.7 ±0.7) x 10 4 
(1.5 ±0.3) x 10 5 
(4.0 ± 1.0) x 10 5 



1.4 ±0.2) x 10 5 
(5.2 ±0.6) x 10 5 
(1.7 ±0.3) x 10 6 
(4.8 ±1.1) x 10 6 
(3.1 ± 1.0) x 10 7 



keep adding one to the power of a randomly chosen variable until the monomial has the 
desired degree. The table show mild exponential growth as degree or dimension grow. This 
was predicted in the theory. 

In conclusion, to integrate a multivariate polynomial, we first algebraically decompose 
each monomial to a sum of powers of linear forms which, as we will see next, can be 
integrated very fast in practice over simplices or over simplicial cones using a few useful 
formulas. Thus we will need a geometric decomposition of our polytopes into those pieces. 

2.2. Valuations and formulas of integration of exponentials over cones and simplices 

We now recall several formulas for the integrals of a power of a linear form over a 
simplex or over a simplicial cone. The idea is that if we can do fast integration for those 
two structures, then we can always rely on two polyhedral decompositions of the input 
polyhedron to obtain the integral. See Subsection 2.5 for details. 

One of the most important properties of integrals over polyhedra is that they can be seen 
as valuations. A valuation F is a linear map from the rational vector space of the indicator 
functions of rational polyhedra P C l d into a rational vector space M. Whenever the 
indicator functions [Pj] of a family of polyhedra Pj satisfy a linear relation £\ [Pj] = 0, 
then the elements P(P«) satisfy the same relation £\ rjP(Pj) = (for a formal definition 
within the polytope algebra, see Chapter 2 of [22]). 

Let C = 5^i=i ^-+ u i be the simplicial cone spanned by linearly independent integral 
vectors Ui,U2, ■■■Ud- The fundamental parallelepiped He of the cone C (with respect to 
the generators u iy i = l,...,d) is the set of points 11c — SiLilA Let us denote by 

vol(IIc) its volume. 

Proposition 1 (Theorem 8.4 in [22]). . There exists a unique valuation I(P)(£) which as- 
sociates to every polyhedron P C V a meromorphic function so that the following properties 
hold 
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(i) If £ is a linear form such that e^ e ' x ' is integrable over P , then 

I(P)(£) = J e {e ' x) dm. 

(ii) For every point s e M. n , one has 

I(s + P)(£) = e^ s) I(P)(£). 

(iii) If P contains a straight line, then I(P) = 0. 

A consequence of the valuation property is the following fundamental theorem. It fol- 
lows from the Brion-Lasserre-Lawrence-Varchenko decomposition theory of a polyhedron 
into the supporting cones at its vertices [19, 22, 21, 20]. 

Lemma 2. Let P be a polyhedron with set of vertices V(P). For each vertex s, let C S (P) 
be the cone of feasible directions at vertex s. Then 

W= E Hs + C s (P))(£). (2) 

Note that the cone C S {P) in Lemma 2 may not be simplicial, but for simplicial cones 
their integrals have explicit rational function formulas. As we see in Proposition 4, one can 
derive an explicit formula for the rational function I(s + C S (P)) in terms of the geometry 
of the cones. 

Lemma 3. Using the valuation property for the valuation I(P)(£) and the linearity over 
the integrands we have that: 

(i) For any triangulation T of the polytope P, we have I(P)(£) = ^ Ae7 - /( A) (/) . 

(ii) For any triangulation T> s of the feasible cone C S (P) at each of the vertices s of the 
polytope P we have I(P)(£) = EseV(P) ^cev s H s + 

Lemma 3 says that if we know how to integrate over simplices or simplicial cones, we 
can integrate over a polytope. We are close to knowing how to do this. By elementary 
integration, and Proposition 1, we have the following. 

Proposition 4. For a simplicial cone C generated by rays u\, 112, ■ ■ ■ Ud (with vertex 0) and 
for any point s 

I(s + C)(£) = vol(n c )e^> n — ^— . (3) 

This identity holds as a meromorphic function of I and pointwise for every £ such that 
(£, Ui) 7^ for all Ui. 
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2. 3. From exponentials to powers of linear forms 

We now consider powers of linear forms instead of exponentials. Similar to I(P), we 
now let L M (P) be the meromorphic extension of the function defined by 

L M (P)(£) = [ (£,x) M dm 



for those £ such that the integral exists. To transfer what we know about integrals of 
exponentials to those of powers of linear forms, we can consider the formula of Proposition 4 
as a function of the auxiliary parameter t: 



[ e<"'*> dm = vol(n c )e<^> f[ — \ 

Js+C i=1 \ t( -i 



—y (4) 



Using the series expansion of the left in the variable t, we wish to recover the value of 
the integral of (£,x) M over the cone. This is the coefficient of t M in the expansion; to 
compute it, we equate it to the Laurent series expansion around t — of the right-hand- 
side expression, which is a meromorphic function of t. Clearly 

1=1 1 ' n=0 1=1 

We say that i is regular if (£,Ui) ^ for every ray of the cone. With this, we can 
conclude the following. 

Corollary 5. For a regular linear form I and a simplicial cone C generated by rays 
Ui, m 2 , . . . u d with vertex s 



L M (s + C)(£) = vol(n c ) ^7\ . (5) 



Ml ((l,*)) M+d 



Otherwise when i is not regular, there is a nearby perturbation which is regular. To obtain 
it, we use £ + e where e = ea is any linear form with a G M n such that (—£ — i, Ui) ^ for 
all Ui, to define a new linear form (depending of a) on the space of meromorphic functions 
in the variable e. Then, applying (5) on the limit as e goes to zero we obtain: 

Since the reader may not be familiar with residues Res e= o(/) and how to calculate them, 
we recall some useful facts on complex analysis (see, e.g., [30] for details). As we observed, 
there is a singularity or pole at e = for a univariate rational function f(e) = ^4 (which 
in this case is explicitly given in Formula (6) of Corollary 5). Recall that if / has Laurent 
expansion f\e) = Yl'kL-m a k ek , the residue is defined as a_i. Given a rational function f(e) 



7 



with a pole at e = there are a variety of well-known techniques to extract the value of the 

p(Q) 

9'(0)' 



residue. For example, if e = is a simple pole (m = 1), then Res e=0 (/) = ^rms- Otherwise, 



when e = is a pole of order m > 1, we can write /(e) = £ m^| £ ) ■ then expand p, q\ in 
powers of e with p(e) = a + ax£ + a 2 e 2 + . . . and qi(e) = b + b\e + b 2 e 2 + . . . . This way 
the Taylor expansion of p(e)/qi(e) at eo is Co + c\E + C2S 2 + C3E 3 + . . . , where Co = f^, and 
Cfe = ^(fl/t — &iCfe_i — b 2 Ck-2 — • • • — &fcCo). Thus we recover the residue Res e= o(/) = c m -\. 
We must stress that the special structure of the rational functions in Corollary 5 can be 
exploited to speed up computation further rather than using this general methodology For 
more on this see [31, 22, 15] and the following discussion. 

Finally, we have all the tools necessary to write down our formula for integration using 
cone decompositions. 

Corollary 6. For any triangulation D s of the tangent cone C S (P) at each of the vertices s 
of the polytope P we have 

L M (p)(i)= J2 lM ( s+c ^- ( ? ) 

2.4- The formula for the simplex 

Suppose now that A C M" is a <i-dimensional simplex (as it may appear in a triangu- 
lation of the polytope P), and i is a linear form on M. n . We say that the linear form £ is 
regular for the simplex A if it is not orthogonal to any of the edges of the simplex. If I is 
regular for A, then it is regular for all tangent cones at each of the vertices. We then find 
the following result special case of Corollary 6. 

Corollary 7 (Brion, see [19]). Let A be a d-simplex with vertices si, . . . , Sd+\ € W 1 . Let 
I be a linear form which is regular w.r.t. A, i.e., (£, Sj) 7^ (£, Sj) for any pair i 7^ j. Then 
we have the following relation. 

r M I / d+1 IP e \ M + d \ 

When £ is regular, Brion's formula is very short; it is a sum of d+l terms. When £ is 
not regular, we can again use a perturbation £ + e where e = ea as in Corollary 5, so that 
the expression of the integral over the simplex reduces to a sum of residues as in (6). 

However, in the special case of a simplex, there exists a computationally more efficient 
method that avoids the calculation of a perturbation a; see [15]. From [15, Theorem 10] 
we find that L AI (A)(£) is the coefficient of the term t M in the Laurent series of the rational 
function 

M! 1 

d vol(A, dm)— ^—r-, (9) 

'(M + dy. n£}(i -t(i, 8j )) 

in the variable i G C. This rational function can be expanded into partial fractions. To 
this end, let K C {1, . . . , d + 1} be an index set of the different poles t — tk '= 1/ (£, Sfe), 
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and fork^K let m k denote the order of the pole, i.e., 

m k = #{ie{l,...,d+l}:(£, = (£, s h ) }. 
Then the rational function can be written as 



Ojfcl Qfc,2 a k,m k 

_| 1_ 



^Al - Sfc ) ^ (1 - t(£, Sk )y ^ ^ (1 - t(£, s fc )) ro * 

where the coefficients are given by certain residues about the pole t = t k . After a 
change of variables, t — t k + e, one obtains the following formula. 

Corollary 8 (Corollary 13 in [15]). Let A be a d- dimensional simplex. Then for an 
arbitrary power (£,x) M of a linear form, we have: 



[ (£, x) M dm = d\ vol(A, dm) . - ^ - V Res £=0 + (t,Sk)) M+d 

To conclude we note that one can even extend the formula above on integrating the 
power of a linear form to the case of a product of powers of several linear forms (see [15]). 

2.5. Should one triangulate or cone decompose? 

One could triangulate the whole polytope, or integrate over each tangent cone. However, 
each cone must be decomposed into simplicial cones. This is the trade-off: we can get away 
with not doing one large polytope triangulation, but we might have to do many smaller 
cone triangulations. 

The number of simplices in a triangulation and the number of simplicial cones con- 
taining in a polytope can significantly differ. Depending on the polytope, choosing the 
right method can determine its practicality. Our experimental results agree with [16] in 
showing that triangulating the polytope is better for polytopes that are "almost simplicial" 
while cone decomposition is faster for simple polytopes. The details will be discussed in 
Section 4. 

Lemma 3 together with the formulas we stated for integration over simplices and cones 
give a general process for computing integrals: 

• We decompose our polynomial as a sum of powers of linear forms. 

• We select a decomposition of the polyhedron in question, either a triangulation or a 
cone decomposition. 

• Apply the formulas to each piece and add up the results via the above results. 
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Figure 1: A pentagon 



2.6. Examples 

2.6.1. Integral values encoded by rational function identities 

Before working out a simple integration example, let us highlight the fact that for 
regular linear forms the integration formulas are given by sums of rational functions which 
we read from the geometry at vertices and possibly a cone decomposition method: Consider 
a pentagon P with vertices (0, 0), (2, 0), (0, 2), (3, 1), and (1, 3) as in Figure 1. 

Then the rational function giving the value of f p (c\x + C2y) M dxdy is 

ns\ /in \M+2 f . nM+2 / . \M+2 / n \M+2 \ 

M! f ( 2 ci) (3CX + C2) (ci + 3c 2 ) (2 c 2 ) \ 

(M + 2)! y Cl (- Cl -c 2 ) ( Cl + c 2 )(2 Cl -2c 2 ) ( Cl + c 2 )(-2 Cl + 2c 2 ) (- Cl -c 2 )c 2 )' 

This rational function expression encodes every integral of the form J p (cix+c 2 ?/) M da; d|/. 
For example, if we let M — 0, then the integral is equal to the area of the pentagon, and 
the rational function simplifies to a number by simple high-school algebra: 

1 ( 4 Cl + 4 ( 3ci + C2 ) 2 + 4 (C1 + 3C2)2 + 4 C2 

2 I -ci - c 2 (ci + c 2 ) (2ci - 2c 2 ) (ci + c 2 ) (-2c! + 2c 2 ) — ci - c 2 

Hence the area is 6. When M and (ci,c 2 ) are given and (ci,c 2 ) is not perpendicular 
to any of the edge directions we can simply plug in numbers to the rational function. For 
instance, when M = 100 and (ci = 3, c 2 = 5) the answer is a fraction with numerator equal 
to 

227276369386899663893588867403220233833167842959382265474194585 
3115019517044815807828554973991981183769557979672803164125396992 

and denominator equal to 1717. When (ci,c 2 ) is perpendicular to an edge direction, we 
encounter (removable) singularities in the rational functions, thus using complex residues 
we can do the evaluation. Note that those linear forms that are perpendicular to some 
edge direction form a measure zero set inside a hyperplane arrangement. 
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(0,1) 


<«>) 


<-i,o) 


,(i,i) 




<l-l\ 




'<0,-l> 








k (0,l) 






<-i,iK ' 














(1,0) 



(0,0) (1,0 



<-i,o> (1,1) 



<0,-l) 




,<o,i> 



(-1,0) (1,0) 



(a) Triangulation method (b) Cone decomposition method 
Figure 2: Example polytopes 



2.6.2. Using the triangulation method 

Take the problem of integrating the polynomial x + y over the triangle A with vertices 
si = (1, 1), s 2 = (0, 1), and s 3 = (1, 0) in Figure 2a. 

The polynomial is already a power of a linear form, and the polytope is a simplex. 
Because £ = (1, 1) is not regular (it is perpendicular to the edge spanned by s 2 and S3), we 
have to build the index set K. Note (£, si) = 2, (£, s 2 ) = 1, and (£, s 3 ) = 1; pick K = {1, 2} 
with mi = 1, m 2 = 2. We now need to compute two values: 

Vertex si. We are not dividing by zero, we can simply plug vectors into Corollary 7, 

(t,3l) 3 



',Sl-S2> 2 

Vertex s 2 : Here, we need to compute a residue. 



Res 



(e + (£,s 2 )) 



1+2 



e=0 



Res 



(e + 1 



,1+2 



e=0 



Finally, j A (x + y) dx dy = 2! x \ x 1! 



si)) 

,,,8-4) = 2/3. 



2.6.3. Using the cone decomposition method 

Next, integrate the polynomial x over the unit square in Figure 2b using the cone 
decomposition algorithm. The polynomial is already a power of a linear form so I = (1, 0). 
The polytope has four vertices that we need to consider, and each tangent cone is already 
simplicial. 

Vertex si = (0, 0): Because {£, si) 1+2 = 0, the integral on this cone is zero. 

Vertex s 2 = (0, 1): For the same reason as Si, the integral on this cone is zero. 

Vertex s 3 = (1,0): At this vertex, the rays are U\ = (0,1), u 2 = (—1,0). Because 
(£, U\) = 0, we need a perturbation vector e so that when i := £ + i, we do not divide by 
zero on any cone (we have to check this cone and the next one). Pick s = (e, e). Then the 
integral on this cone is 
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M' (l + e) 1+2 V 

* vol(n c ) Res e=0 ; \,! r = x 1 x -2 = -2/6. 



(M + d)\ v w £ " u e (- e )(l + e) (1 + 2)! 

Vertex s 4 = (1, 1): The rays are U\ = (— 1,0), m 2 = (0,-1). Again, we divide by zero, 
so we perturbate i by the same e. The integral on this cone is 

M' fl + 2rl 1+2 1' 

' vol(n c ) Res £=0 ,T, c = t — — tt x 1 x 5 = 5/6. 



(M + d)! v w e(e)(l + e) (1 + 2)! 

The integral f p x dx dy = + — 2/6 + 5/6 = 1/2 as it should be. 



3. How the software works 

LattE was originally developed in 2001 as software to study lattice points of convex 
polytopes [32] . The algorithms used combinations of geometric and symbolic computation. 
Two key data structures are rational generating functions and cone decompositions, and 
it was the first ever implementation of Barvinok's algorithm. LattE was improved in 2007 
with various software and theoretical modifications, which increased speed dramatically. 
This version was released under the name LattE macchiato; see [33]. Now in 2011, our new 
release LattE integrale has extended its capabilities to include the computation of exact 
integrals of polynomial functions over convex polyhedra. The new integration functions 
are C++ implementations of the algorithms provided in [15] with additional technical 
improvements (including an important new set of data structures for the manipulation of 
truncated series). A key distinction between LattE integrale and other software tools 
is that our algorithms give the exact evaluation of the integral since our implementation 
uses exact rational arithmetic. The code of this software is freely available at [14] 

The new implementation of LattE integrale allows us to calculate the integral of a 
sum of powers of linear forms over an arbitrary polytope. Alternatively, we can calculate 
the integral of a sum of monomials by decomposing each monomial into a sum of powers 
of linear forms using Formula (1), then integrating these powers of linear forms. 

This section starts with a discussion of our new data structure for manipulating poly- 
nomials and linear forms, then we describe the format LattE integrale expects for the 
input polytopes, and we end with a detailed explanation of the two main algorithms. 

3.1. Input format and data structures 

The input format for the polynomials is identical to that of the Maple programs of [15]: 

• A polynomial is represented as a list of its monomials in the form 

[monornialx , monomial 2 ,...], 

where monorniali is represented by 

[coefficient , [exponent-vector] ] . 
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For example, 3xgxf x\ + 7x\ x\ is input as [ [3 , [2 , 4 , 6] ] , [7 , [0 , 3 , 5] ] ] . 
• To deal directly with powers of linear forms, the input format is 

[lineai — terrni, lineai — term 2, ■ . .], 
where lineai — terrrii is represented by 

[coefficient , [power, [coefficient-vector]]] . 

For example, 3(2x + 4a* + 6a;2) 10 + 7(3a;i + 5a; 2 ) 12 is input as [[3, [10, [2,4,6] ]] , 
[7, [12, [0,3,5]]]]. 

In [15], the integration over simplices was first implemented in Maple, and so there 
was no control over the data structures used to store the data. We have implemented 
the simplex integration algorithm in C++ with a sophisticated data structure and have 
developed a new algorithm that integrates over the tangent cones of a polytope. Currently, 
we are using burst tries, a data structure designed to have cache-efficient storage and search, 
due to the fact that they are prefix trees with sorted arrays of stored elements as leaves [34]. 
Such a data structure is performance-critical when computing residues, as a comparison 
with a linked-list implementation showed. In our implementation, each node corresponds 
to a particular dimension or variable and contains the maximal and minimal values of the 
exponent on this dimension. The node either points to another node a level deeper in the 
tree or a list of sorted elements. 

The input rational polytope P could be given to LattE integrale by an h-representation 
or v-representation. The input format is the same as in previous versions and it is ex- 
plained in [14]. Although the theory we presented earlier works for both full- dimensional 
and non-full dimensional rational polytopes, the current release of LattE integrale is 
only guaranteed to do integration and volume computation in full-dimensional polyhedra. 
It is worth stressing the old capabilities for counting lattice points still work for polytopes 
of all dimensions and we impose no arbitrary limit on the size or dimension of the input. 
LattE integrale relies on Cddlib [35] or 4ti2 [36] for all basic polyhedral calculations 
such as computation of dimension. 

Our data structures are specialized for polytopes with vertices of integer coordinates. 
In order to integrate over rational polytopes, we first dilate them and perform a change of 
variables. If P is a <i-dimensional rational polytope and aP is a dilation by a > that 
makes P integer, then our software operates on the vertices of aP and rescales the final 
integral by the following well-known change of variables: 



After this transformation, we apply Formula (1) to transform the polynomial into powers 
of linear forms. 

When integrating polytopes other than simplices, there are two options based on the 
formulas presented in Section 2: (i) Triangulate the polytope and apply the algorithm for 
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each simplices individually, or (ii) Triangulate each tangent cone and integrate each one 
using the cone decomposition algorithm. Therefore, a key step in all our computations is 
to find either a triangulation of the polytope or a triangulation of each of its tangent cones. 
Once more, this step relies on Cddlib or 4ti2, because when we triangulate we compute 
a regular triangulation via a convex hull [32, 33]. We now explore the two integration 
algorithms in greater detail. 

3.2. Integrating powers of linear forms by polytope triangulation 

After we decompose the polynomial to a sum of powers of linear forms and after finding 
a triangulation of the polytope, Algorithm 1 loops over these two sets and integrates each 
linear form/simplex pair individually using Corollaries 7 and 8. 

Algorithm 1 Integrate using polytope triangulation 

Input: F = c j(^ji x) Mj , sum of powers of linear forms. 
Input: P, a full-dimensional polytope. 
Output: integral of the linear forms F over the polytope P. 
integral <— { integral is a rational data type} 
let T be a list of simplices that form a triangulation of P 
for all simplices A in T do 

for all linear forms c(£,x) M in F do 
if £ is regular on A then 

integral integral + c x integrateSimplexRegular(£, M, A) 
else 

integral t— integral + c x integrateSimplexResidue(£, M, A) 
end if 
end for 
end for 

return integral 



In Algorithm 1, the linear forms are represented as a burst trie, the triangulations are 
stored in a linked list, and each simplex is a simple two-dimensional array containing the 
vertices Sj, . . . , sj+i- 

When £ is regular on A, the integrateSimplexRegular function plugs in numbers and 
vectors into Corollary 7. Also, the terms in the numerator and denominator are evaluated 
in a rational data type, and so no floating-point divisions are performed. 

When £ is not regular, the integrateSimplexResidue function (Algorithm 2) applies 
Corollary 8 and must find an index set K C {1, . . . , d + 1} of different poles t = l/(£, 
and compute \K\ residues. Let k 6 K and let denote the order of the pole, i.e., 

m k = #{ie{l,...,d+l}:(£,8i) = (£, s k ) }. 

The problem has now been reduced to evaluating Formula (11). 
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£ m k Yl ( e + (£, s k - Si))"* L ^ 1 n(^+ (t^k-Si))™*' 

i£K ieK 

where [e a ]p means the coefficient of e a in the Laurent series of expression p. 

To compute Formula (11), we expand the polynomial in terms of e in the numerator 
truncated to degree m k — 1. We then find the first m k terms in the polynomial expansion 
of each l/(e + (£, s k — Si)) m \ i ^ k term using the general binomial theorem. To make the 
notation easy, let b = (£, s k — Sj) G Z, then the degree m k — 1 polynomial of (e + b)~ mk is 

p(e) = a e%- m * + a^b^' 1 + ■■■ + a^^e^- 1 ^"" 1 ^ 1 , a 3 - = ( m * +J ~ 1 ] (_i)i. 

This is a polynomial in e with rational coefficients. For efficiency reasons, we factor 
p(e) = bm .+ mk _i_ p(e), p G Z[e]. The integrateSimplexResidue and truncatedMultiply func- 
tions both implement these ideas. 

Algorithm 2 integrateSimplexResidue 

Input: coefficients of a linear form and M, integer power. 

Input: A, simplex with vertices Si, . . . , Sd+i- 

Output: The integral of (£,x) M over A. 

Let pi <— 1,P2 <— 1 be polynomials in e. 

Let r/ <- 1 be a rational data type. 

Let sumResidue be a rational data type. 

Make the index set K of unique poles. 

for all k in K do 
rf <- 1 

Pi the expansion of (e + (£, Sk)) M+d up to degree — 1 {pi G Z[e]} 
for all i in K and i 7^ k do 

rf ^ rf x (£, Sk )-(m+m k -i) 

p 2 -f- the expansion of (e+ (£, s k — Si))~ mi up to degree m k — 1 with (£, s k )~( mi+mk ~ 1 > 
factored out. {p 2 G Z[e]} 
Pi ^— truncatedMultiply (pi, p2,rrik — 1) 
end for 

Let c be the coefficient of the degree m k — 1 term in pi(e). 
sumResidue sumResidue + rf x c 
end for 

return abs(det(s\ — sa+i, ■ ■ ■ , — Sd+i)) x n^g)i x sumResidue 



Finally, Algorithm truncatedMultiply (p,q,k) takes two polynomials p, g G Z[e] and re- 
turns their product up to and including terms of degree k. Our implementation is very 
simple (e.g., not using any special multiplication algorithms) but the cache-efficient use of 
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the burst tries leads to speed-ups wheatn compared to a naive implementation with arrays. 
We note that asymptotically faster multiplication algorithms exists (see, e.g., [37]), which 
might lead to further improvements. 

3. 3. Integrating powers of linear forms by cone decomposition 

After triangulating each tangent cone into simplicial cones, the computation is very 
similar to the polytope-triangulation case: if I is regular on the rays of the cone, we plug in 
values into Corollary 5, else we perturb £ and perform a residue calculation. Algorithm 3 
implements this idea. 

Algorithm 3 Integrate using the cone decomposition method 

Input: F = c j (^ji x ) Mj > powers of linear forms. 
Input: P, a full-dimensional polytope. 
Output: integral of the linear forms over P. 

integral <— { integral is a rational data type} 

Let C be a list of tangent cones P. 

Make T be a list of triangulated cones in C. {A cone in T is in the form (s; ui, . . . , ua), 
where s is a vertex and Ui are rays} 
for all linear forms c(£,x) M in F do 

Let R C T be cones that £ is regular on. 

for all (s; Mi, . . . , Ud) in R do 

integral -t- integral + c x integrateConeRegular(£, M, s,u\, . . . , Ud) 

end for 

Pick e = e{a\, . . . , a<j) where ^ 6 Z so that we do not divide by zero, 
for all (s; u\, . . . , u d ) in T \ R do 

integral integral + c x integrateConeResidue(£, M, e, s,ui, . . . , Ud) 
end for 
end for 

return integral 



In integrateConeResidue, £ is perturbed by setting £ := £+i, where e is a vector in terms 
of e with coefficients picked on the moment curve with alternating signs. We repeatedly 
pick a random t G Z + and set ii = t 1 ^ 1 {—l) 1 ^ 1 e for i = 1,2, ... ,d until (—(£ + i),Ui) 
is non-zero for every simple cone at every vertex in Corollary 5. Then the residues are 
computed using the general binomial theorem and truncated series multiplication like in 
integrateSimplexResidue. 

3.4- A special case: computing volumes 

Computing the volume of a polytope is equivalent to integrating the monomial 1 over 
the polytope. We again have the two same options when computing volumes as we did 
when computing integrals. Instead of using the algorithms above, we can simplify the 
computation. In the triangulation based approach, we find a triangulation of the polytope 
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and sum the volume of each simplex. The volume of a specific simplex is obtained by 
taking a determinant. In the cone decomposition approach, we triangulate each tangent 
cone and apply Corollary 5 with M = and any random vector I. If we do divide by zero, 
instead of finding residues, we simply pick a new random I and start the computation over. 

4. Experiments 

We did thorough testing of the implementation and performed new computational 
benchmarks. We report on four different test classes: 

1. We expand the computational limits for integrating over simplices described in [15]. 

2. Next, we integrate random monomials over three families of polytopes: (1) simple, 
(2) simplicial, and (3) neither simple nor simplicial. 

3. Because our volume methods are optimized versions of the integration methods, we 
also compute the volumes of the same polytopes in the last case above. 

4. Finally, we compare LattE integrale to other software tools and computational 
studies [16, 38, 39, 3]. 

Our integration and volume experiments input data, along with running times and the 
results of integration and volumes are available on the LattE website [14]. 

4-1. Integration over simplices 

In [15], the theory of integration over simplices was developed and a fair amount of 
Maple experiments were carried out to show the potential of the methods. In this section, 
the experiments we performed clearly indicate that this C++ implementation is at least 
two orders of magnitude faster than the preceding Maple code; compare Tables 5 and 6 in 
[15] with Tables 2 and 3 in this paper. In Table 2, we used Maple to generate powers of 
random linear forms and randomly generated simplices. The coefficients of each linear form 
were picked uniformly over [0, 100] D Z. We did the integration using LattE integrale. 

Next, in Table 3, we used Maple to generate monomials and simplices. We again 
integrate using LattE integrale. We measure time from the start of program execution 
to termination, which includes file I/O, system calls, child process time, the time to find 
tangent-cones, and triangulation time. All triangulations were computed with the software 
package cddlib version 0.94f [35]. All computations were performed on a 64-bit Ubuntu 
machine with 64 GB of RAM and eight Dual Core AMD Opteron 880 processors. We 
applied a 600-second maximum running time to this program; tasks taking longer are not 
benchmarked. 
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Table 2: Average integration time plus or minus one standard deviation when integrating 
one power of a linear form over a random rf-simplex (in seconds over 50 random forms) 



Exponent M 



d 


2 


10 


20 


50 


100 


300 


1000 


10 


0.00 ±0.00 


0.00 ±0.00 


0.00 ±0.00 


0.00 ±0.00 


0.00 ±0.00 


0.00 ±0.00 


0.01 ±0.00 


20 


0.00 ±0.00 


0.00 ±0.00 


0.00 ±0.00 


0.00 ±0.00 


0.00 ±0.00 


0.00 ±0.00 


0.01 ±0.00 


50 


0.00 ±0.00 


0.00 ±0.00 


0.00 ±0.00 


0.00 ±0.00 


0.00 ±0.01 


0.01 ±0.00 


0.03 ±0.01 


100 


0.02 ±0.00 


0.02 ±0.00 


0.02 ±0.00 


0.02 ±0.00 


0.02 ±0.00 


0.03 ±0.00 


0.09 ±0.00 


300 


0.35 ±0.01 


0.36 ±0.01 


0.36 ±0.01 


0.36 ±0.01 


0.38 ±0.01 


0.42 ±0.01 


0.66 ±0.02 


400 


0.78 ±0.02 


0.79 ±0.03 


0.79 ±0.02 


0.80 ±0.03 


0.82 ±0.03 


0.90 ±0.03 


1.25 ±0.04 



Table 3: Average and standard deviation of integration time in seconds of a random 
monomial of prescribed degree by decomposition into linear forms over a d-simplex (average 
over 50 random forms) 



Degree 



a 


1 
1 


o 
z 





i n 
1U 


on 


oU 


/in 


OU 


1 nn 
1UU 


zUU 


oUU 


2 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


0.1 


1.0 


3.8 




0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


0.1 


0.4 


1.7 


3 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


0.1 


0.2 


2.3 


38.7 


162.0 




0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 


0.1 


1.4 


24.2 


130.7 


4 


0.0 


0.0 


0.0 


0.0 


0.0 


0.1 


0.4 


0.7 


22.1 








0.0 


0.0 


0.0 


0.0 


0.0 


0.1 


0.3 


0.7 


16.7 






5 


0.0 


0.0 


0.0 


0.0 


0.1 


0.3 


1.6 


4.4 










0.0 


0.0 


0.0 


0.0 


0.0 


0.2 


1.3 


3.5 








6 


0.0 


0.0 


0.0 


0.0 


0.1 


1.1 


4.7 


15.6 










0.0 


0.0 


0.0 


0.0 


0.1 


1.0 


4.3 


14.2 








7 


0.0 


0.0 


0.0 


0.0 


0.2 


2.2 


12.3 


63.2 










0.0 


0.0 


0.0 


0.0 


0.2 


1.7 


12.6 


66.9 








8 


0.0 


0.0 


0.0 


0.0 


0.4 


4.2 


30.6 


141.4 










0.0 


0.0 


0.0 


0.0 


0.3 


3.0 


31.8 


127.6 








10 


0.0 


0.0 


0.0 


0.0 


1.3 


19.6 














0.0 


0.0 


0.0 


0.0 


1.4 


19.4 












15 


0.0 


0.0 


0.0 


0.1 


5.7 
















0.0 


0.0 


0.0 


0.0 


3.6 














20 


0.0 


0.0 


0.0 


0.2 


23.3 
















0.0 


0.0 


0.0 


1.3 


164.8 














30 


0.0 


0.0 


0.0 


0.6 


110.2 
















0.0 


0.0 


0.1 


4.0 


779.1 














40 


0.0 


0.0 


0.0 


1.0 


















0.0 


0.0 


0.3 


7.0 
















50 


0.0 


0.0 


0.1 


1.8 


















0.0 


0.1 


0.5 


12.9 
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4-2. Integration over general polytopes 

We tested the triangulation and cone decomposition integration methods on polytopes 
and their duals across dimension, vertex counts, and over monomials with different degrees. 
For each polytope dimension and vertex count we constructed 50 random polytopes by 
taking the convex hull of random points using Polymake [40]. For primal polytopes of 
dimension d, the number of vertices considered goes from d + 2 to d + 25. When zero 
is not in the interior of the polytope, we translated the centroid to the origin before 
constructing the dual polytope. Then we integrated each polytope and its corresponding 
dual polytope over a new random monomial of a set degree. Because of the construction 
method, most primal polytopes are simplicial and the duals are mostly simple polytopes. 
We also integrated over G. M. Ziegler's database of polytopes [41], which contains polytopes 
that are not simplicial nor simple. 

Simple and simplicial polytopes We present the results for dimensions 5, 6, and 7. We 
tested both algorithms on the primal polytopes starting with their v-representation. 
For their duals we tested the triangulation and cone decomposition methods starting 
from their h-representations. We also did experiments in dimension 3 and 4 but the 
numbers are too close to each other to show a clear trend of which is the fastest 
method (see [42]). We only report those test classes for which every polytope in 
the test class finished under 600 seconds for both the triangulation and the cone- 
decomposition method. 

In Figures 3, 4, and 5, we display histograms that on the x-axis plot the time difference 
between the two integration methods, the |/-axis shows the degree of monomials and 
the z-axis presents the number of random polytopes (in the respective dimensions 
5, 6, 7). A particular solid bar in position (a*,b*) tallies the number of random 
polytopes for which the time difference between the two algorithms was a* seconds 
when integrating a a monomial of degree b*. We define the time difference as the time 
taken by the triangulation method minus the time taken by the cone-decomposition 
method. There is one histogram for the primal polytopes and one for the dual 
polytopes. For example, in Figure 3 on the top (primal), there are eight colors on the 
bars, one for each degree (corresponds to a row). The bars in the degree 50 row, show 
the average time differences are always less than or equal to zero. This shows that 
the cone method was slower in the majority of simplicial problems. For comparison, 
on the bottom (dual), most of the mass has positive time difference, which indicates 
the cone method wins over the integration of higher degree monomials. More tables 
will be available online [42]. 

In conclusion, our results for dimension higher than four show integrating monomials 
coincide qualitatively with the observation of [16] for volume computation (polyno- 
mial of degree zero): the triangulation method is faster for simplicial polytopes (mass 
on histograms is concentrated on negative time differences) while the cone decompo- 
sition is faster for simple polytopes (mass on histograms is concentrated on positive 
time differences). The trends are very clear while in dimension less than four, the 
timings are too close to each other to give a clear-cut trend. 
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Zero-one polytopes In Figure 6, we present another histogram comparing the cone de- 
composition and triangulation methods on Ziegler's database of polytopes [41], which 
contains many zero-one polytopes and a few other polytopes. We translate each poly- 
tope so that its centroid is the origin, thus its dual is well-defined. Then for each 
polytope and its dual, we integrate 50 random polynomials of a set degree. We 
skipped non-full-dimensional polytopes and a few others that did not finish within 
30 minutes. The figure displays the histogram of the differences in running times 
between the two integration methods for monomials of eight different degrees (1, 2, 
5, 10, 20, 30, 40, and 50). As before, the z-axis is frequency, x-axis is the range of 
the time differences, and y-axis gives the degree used. For example, there are over 
200 cases where integrating a degree 40 monomial over the dual polytopes took the 
triangulation method 90 seconds more than the cone-decomposition method. 

The behavior we observed before for simple vs. simplicial polytopes still mostly holds 
for these tests, except we see that the two methods finish within 20 seconds of each 
other most of the time (mass of histogram is centered at 0). The variation is then 
not as strong as before. 

4-3. Volume Experiments 

Volume computation is an important special case of integration that has received atten- 
tion by several researchers, thus we also tested the triangulation and cone decomposition 
methods on the same database of random polytopes and their duals, and on Ziegler's 
database to see the performance of volume evaluation. 

Simple and simplicial polytopes Each test class contains 50 polytopes for each dimen- 
sion and we only considered tests where both methods finished within 600 seconds 
for the same polytope. While the triangulation method is still faster for simplicial 
polytopes and the cone-decomposition method is faster for simple polytopes, the his- 
tograms in Figure 7 suggest that both methods finish quite close to each other in 
small dimension. When dimension starts growing there is a more pronounced dif- 
ference between the methods (i.e., the mass of the histogram is more spread toward 
positive or negative values of the time difference). 

Zero-one polytopes In Table 4 and 5, we apply the triangulation and cone decomposi- 
tion volume methods to Ziegler's database [41] and their duals. If a polytope did not 
contain the origin, we centered it so that its dual is defined. Again, we skipped non- 
full-dimensional polytopes and a few others that did not finish within 30 minutes. 
Faster timings are shown in bold. When computing volumes of primal polytopes in 
Ziegler's database, triangulation is faster more often. But for finding the volume of 
the dual polytopes there is no clear faster method. 
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Dimension 5: Primal Polytopes (mostly simplicial) 
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Dimension 5: Dual Polytopes (mostly simple) 




U r£7 




ungaaQQQnQnn 



/ dflgrss 40 
J degrsfl 30 
/ dflgrofl 20 
/ d«gr« 10 
/ degrsfl 5 
degree 2 



3 10 13 SO S3 30 33 40 43 30 33 00 03 70 73 SO S3 80 B3 100 

TirfB Difference (Cone decomposiion Faster) — > 



V degree I 



Figure 3: Histogram of the time difference between the triangulation and cone- 
decomposition methods for integrating over random polytopes in dimension 5 
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Dimension 6: Primal Polytopes (mostly simplicial) 




JS 40 

30 



-20 - I B - 1 5 - 1 7 - 1 6 - I E - I 4 - 1 3 - I fi - I I -10 -5 -3 -7 -0 -5 -4 -3 -fi - I 



— (Triangula ion Faster) Tine DfFe rents 



Dimension 6: Dual Polytopes (mostly simple) 
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Figure 4: Histogram of the time difference between the triangulation and cone- 
decomposition methods for integrating over random polytopes in dimension 6 
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Dimension 7: Primal Polytopes (mostly simplicial) 



E 

u 

3- 



100 




an n a n a QQ3 q &a a & 

nDQB 
mQ 



100 -55 -50 -35 -30 -75 -70 -05 -00 -55 -50 -45 -40 -35 -30 -25 -fiO - 15 

< — (Trisn^jIsiionFsster) Time Difference 




Dimension 7: Dual Polytopes (mostly simple) 
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Figure 5: Histogram of the time difference between the triangulation and cone- 
decomposition methods for integrating over random polytopes in dimension 7 
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Ziegler's Database: Primal Polytopes 
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Ziegler's Database: Dual Polytopes 
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Figure 6: Histogram of the time difference between the triangulation and cone- 
decomposition methods for integrating over the polytopes in Ziegler's Database 
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Volume: Primal Polytopes (mostly simplicial) 
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Figure 7: Histogram of the time difference between the triangulation and cone- 
decomposition methods for finding the volume of random polytopes 
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Table 4: The triangulation vs. cone-decomposition method for finding volumes in Ziegler's 
Database: Part I 









Primal 






Dual 










Time 


(sec) 




Time (i 


sec) 


Polytope 


Dim. Vertices 


Cone. 


Triang. 


Vertices 


Cone. 


Triang. 


3simp3simp.vrep.latte 


6 


44 


5.61 


6.10 


32 


1.11 


1.15 


cyclic_4_8.vrep.latte 


4 


8 


0.09 


0.06 


20 


0.02 


0.10 


neighborly _4_8.vrep.latte 


4 


8 


0.12 


0.03 


20 


0.03 


0.06 


SharirCube.vrep.latte 


3 


8 


0.03 


0.03 


6 


0.11 


0.02 


HC_6-32.vrep.latte 


6 


32 


2.29 


2.06 


44 


3.25 


3.22 


HC_7-64.vrep.latte 


7 


64 


13.42 


75.85 


78 


61.68 


762.12 


HC_8-128.vrep4atte 


8 


128 


85.85 




144 


15007.50 




MJ_16-17.vrep4atte 


16 


17 


2.60 


2.48 


17 


0.07 


0.04 


OA_5-10.vrep.latte 


5 


10 


0.22 


0.08 


22 


0.18 


0.11 


OA_5-18.vrep4atte 


5 


18 


0.46 


0.32 


19 


0.34 


0.10 


OA_5-24.vrep.latte 


5 


24 


0.81 


0.58 


18 


0.22 


0.13 


OA_6-13.vrep4atte 


6 


13 


0.53 


0.20 


56 


0.52 


5.37 


OA_7-18.vrep.latte 


7 


18 


3.36 


0.82 


146 


13.96 


1827.83 


OA_8-25.vrep.latte 


8 


25 


38.55 


10.44 


524 


4116.93 




OA_9-33.vrep.latte 


9 


33 




648.77 


1870 






AS_6-18.vrep4atte 


6 


18 


1.26 


0.52 


121 


1.40 


65.63 


BIR3_4-6.vrep4atte 


4 


6 


0.12 


0.02 


9 


0.01 


0.01 


BIR4_9-24.vrep4atte 


9 


24 


6.26 


2.22 


16 


1.42 


0.17 


BIR5_16-120.vrep4atte 


16 


120 






25 




488.78 
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Table 5: The triangulation vs. cone-decomposition method for finding volumes in Ziegler's 
Database: Part II 









Primal 






Dual 










Time (; 


sec) 




lime 


(sec) 


Polytope 


Dim. Vertices 


Cone. Triang. 


Vertices 


Cone. 


Triang. 


CF_10-ll.vrep.latte 


10 


11 


0.37 


0.33 


11 


0.03 


0.00 


CF_4-5.vrep.latte 


4 


5 


0.01 


0.02 


5 


0.02 


0.02 


CF_5-6.vrep.latte 


5 


6 


0.04 


0.03 


6 


0.00 


0.00 


CF_6-7.vrep.latte 


6 


7 


0.04 


0.05 


7 


0.01 


0.00 


CF_7-8.vrep.latte 


7 


8 


0.08 


0.09 


8 


0.02 


0.02 


CF_8-9.vrep.latte 


8 


9 


0.14 


0.13 


9 


0.02 


0.01 


CF_9-10.vrep.latte 


9 


10 


0.22 


0.20 


10 


0.01 


0.02 


CRO_3-6.vrep.latte 


3 


6 


0.04 


0.01 


8 


0.00 


0.02 


CRO_4-8.vrep.latte 


4 


8 


0.12 


0.05 


16 


0.00 


0.03 


CRO_5-10.vrep.latte 


5 


10 


0.17 


0.10 


32 


0.01 


0.33 


CUT3_3-4.vrep.latte 


3 


4 


0.00 


0.01 


4 


0.01 


0.01 


CUT4_6-8.vrep.latte 


6 


8 


0.16 


0.06 


16 


0.00 


0.05 


CUT5_10-16.vrep.latte 


10 


16 


2.72 


0.94 


56 


38.38 


2046.74 


CYC_5-8.vrep.latte 


5 


8 


0.11 


0.05 


20 


0.00 


0.10 


EG_5-12.vrep.latte 


5 


12 


0.32 


0.13 


40 


0.15 


0.79 


EQU_5-7a.vrep.latte 


5 


7 


0.06 


0.05 


10 


0.01 


0.04 


EQU_5-7b.vrep.latte 


5 


7 


0.09 


0.05 


10 


0.01 


0.00 


HAM_8-16.vrep.latte 


8 


16 


1.57 


0.60 


256 


0.15 




HC_3-4.vrep.latte 


3 


4 


0.03 


0.00 


4 








HC_4-8.vrep.latte 


4 


8 


0.15 


0.05 


16 


0.02 


0.04 


HC_5-16.vrep.latte 


5 


16 


0.48 


0.24 


26 


0.32 


0.25 


CNG_5-6a.vrep.latte 


5 


6 


0.04 


0.03 


6 


0.02 


0.01 


MJ_32-33.vrep.latte 


32 


33 


82.90 


83.86 


33 


1.44 


0.14 


CNG_5-6b.vrep.latte 


5 


6 


0.02 


0.03 


6 


0.01 


0.02 
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4-4- Comparison to other software 

There are two general classes of algorithms for finding volumes and integrals over poly- 
topes: numerical and exact. Numerical algorithms approximate the valuation on the poly- 
tope and involve error bounds, whereas exact algorithms do not contain a theoretical error 
term. However, exact algorithms may contain errors when they use finite digit integers 
or use floating-point arithmetic. In order to sidestep this problem, LattE integrale uses 
NTL's arbitrary length integer and rational arithmetic [43] compiled with the GNU Multi- 
ple Precision Arithmetic Library [44]. The obvious downside to exact arithmetic is speed, 
but this cost is necessary to obtain exact answers. In this section, we compare our ex- 
act algorithms with other software tools and algorithms that use numerical algorithms or 
non-exact arithmetic. 

Vinci contains different algorithms for finding polytope volumes and in fact imple- 
mented the same decompositions we used in our software (see [16]). We tested against 
Vinci 1.0.5, and Table 6 compares LattE's cone decomposition method with Vinci's fastest 
method HOT (Hybrid Orthonormalisation Technique). We ran LattE's cone decomposi- 
tion method starting from the h-representation. Because the HOT method requires both 
an h- and v-representation of the polytope, we also report the time used by CDD [26] to 
convert an h-representation to a v-representation. We also break down time spent in LattE 
for finding the vertices, finding the rays at each vertex, triangulation, and the time spent 
in the main cone decomposition integration method. 

It is clear that the HOT method is faster and usually accurate when applied on the 
Vinci database (these polytopes are available from [16]), but because of non-exact arith- 
metic, it can give incorrect results. In fact we found that Vinci's cone decomposition 
method contained a bug: Vinci's cone decomposition method found the correct volumes 
for the cubes and random-hyperplane polytopes, but reported incorrect or negative vol- 
umes for most polytopes in the Vinci database. We also explored how well Vinci can 
compute volumes of polytopes where each vertex contains small and large positive num- 
bers. In Table 7, we tested the accuracy of Vinci's HOT method on cyclic polytopes. We 
constructed these d- dimensional polytopes by taking the the convex hull of k + d points 
(t, t 2 , t 3 , . . . , t d ) G Z d for t = 5, 6, . . . , 5 + k + d — 1. For very small dimensions, the HOT 
method does well, but gives incorrect or zero volumes already in dimension six. 

Another comparison we made was to the paper [3], where it is claimed that exact vol- 
umes are computed by integration. The authors report seven volumes for different poly- 
topes. LattE integrale's triangulation and cone decomposition method agrees with their 
calculations except in the last case. For P-j the correct volume is 1 /622080 ~ 1. 607510 xlO -6 
but they calculate 1.56439 x 10 -6 . Presumably, because of non-exact arithmetic, their an- 
swer has only one digit of accuracy. 

4-5. Numerical methods 

M. Korenblit and E. Shmerling present a numerical integration algorithm in [39] which 
is based on a special decomposition of the integral into regions that have well-defined upper 
and lower limits of integration that, on an ordering of the variables, xi,£2, . . . ,Xd, x% is 
expressed only in terms of X\, . . . , It is known that achieving such a decomposition 
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Table 6: Time breakdown between LattE integrale's cone-decomposition and Vinci's 
HOT method with CDD 



Vinci LattE 



Polytope 


HOT 


Cddlib 


Vertices 


Rays 


Triang. 


Cone. 


cube-9 


0.03 


0.08 


0.02 1 


0.06 


0.02 


0.02 


cube- 10 


0.11 


0.18 


0.04 1 


0.15 


0.02 


0.06 


cube-14 


141.65 


7.99 


1.24 1 


4.67 


0.69 


1.26 


rh-8-20 


0.13 


0.89 


0.11 1 


0.49 


0.04 


7.00 


rh-8-25 


0.43 


2.63 


0.32 1 


1.14 


0.14 


80.25 


rh-10-20 


0.96 


2.21 


0.25 1 


1.80 


0.14 


98.26 


rh-10-25 


5.71 


12.49 


1.07 1 


8.80 


0.44 


3989.25 


CC 8 (9) 


0.04 


0.22 


0.07 1 


0.12 


0.39 


0.40 


CC 8 (10) 


0.08 


0.52 


0.16 1 


0.22 


0.97 


0.88 


CC 8 (11) 


0.18 


1.18 


0.03 1 


0.42 


1.84 


1.76 


ccp 5 


0.00 


0.07 


0.09 2 


0.00 


0.10 


0.09 


cross 8 


0.00 


0.39 


0.50 2 


0.00 


0.06 


0.04 


cross 9 


0.00 


1.57 


2.15 2 


0.00 


0.12 


0.11 


rv-8-10 


0.00 


0.08 


O.OO 1 


0.03 


0.02 


0.00 


rv-8-11 


0.00 


1.99 


0.08 1 


0.23 


0.03 


0.01 


rv-10-12 


0.00 


0.12 


0.02 1 


0.09 


0.04 


0.01 


rv-10-14 


0.00 


1061.49 


29.50 1 


64.96 


0.10 


0.07 



1 Computed with 4ti2. 

2 Computed with Cddlib. 
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Table 7: Comparison between LattE integrale and Vinci on finding the volume of cyclic 
polytopes 









k 


Dim. Tool 1 


2 


3 4 5 


2 


LattE 1 
Vinci 1 


4 
4 


10 20 35 
10 20 35 


3 


LattE 2 

Vinci 1.999999999988 


16 70 224 588 
15.99999999999 69.99999999991 224.0000000006 587.9999999986 


4 


LattE 12 

Vinci 11.99999993201 


192 1512 8064 33264 
191.9999999913 1511.99999999 8063.999999892 33263.99999989 


5 


LattE 288 9216 133056 1216512 8154432 
Vinci 287.9996545868 9216.000252236 133055.9883262 1216511.998301 8154431.872519 


6 


LattE 34560 
Vinci 34561.951223 


2211840 59304960 948879360 10600761600 
1935359.822684 58060819.63341 885910920.3761 10336274212.34 


7 


LattE 24883200 
Vinci 25744201.0524 


3185049600 



160123392000 4554620928000 86502214656000 




is equivalent to the so-called Fourier-Motzkin elimination procedure [45] and as such it is 
of exponential complexity. The paper [39] gives an application to finding the probability a 
random-coefficient polynomial has one or two real roots in the interval [—1, 1]. To do this, 
they use their software to find the volume of a polytope. They calculate 2.79167; however, 
we verified that the correct volume is 31/12 = 2.583 which gives their method one digit of 
accuracy. 

A more interesting comparison is to CUBPACK, a Fortran 90 library which estimates the 
integral of a function (or vector of functions) over a collection of rf-dimensional hyper- 
rectangles and simplices [38]. This comparison is very interesting because CUBPACK uses an 
adaptive grid to seek better performance and accuracy. All integration tests with CUBPACK 
in dimension d were done with a product of linear forms with a constant term over a random 
(i-dimensional simplex where the absolute value of any coordinate in any vertex does not 
exceed 10. For example, we integrated a product of inhomogeneous linear forms such as 
(l + 2x-j^y)(2-5x) over the simplex with vertices (10, 0), (9, 9), (1, 1). In Table 8, LattE 
was run 100 times to get the average running time, while CUBPACK was run 1000 times due 
to variance. Both the dimension and number of linear forms multiplied to construct the 
integrand were varied. 

As shown in Table 8, LattE integrale tends to take less time, especially when the 
number of forms and dimension increases. The table does not show the high variance that 
CUBPACK has in its run times. For example, the 5-dimensional test case with 6 linear forms 
had a maximum running time of 2874.48 seconds, while the minimum running time was 0.05 
seconds on a different random simplex. This contrasted starkly with LattE integrale, 
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Table 8: Average Time for LattE integrale and CUBPACK for integrating products of 
inhomogeneous linear forms over simplices. 













Number of linear factors 








d 


Tool 


1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


2 


LattE 
CUBPACK 


0.0001 

0.0027 


0.0002 

0.0014 


0.0005 

0.0016 


0.0008 

0.0022 


0.0009 

0.0064 


0.0019 

0.0052 


0.0038 
0.0014 


0.0048 
0.0002 


0.0058 
0.0026 


0.0089 

0.0213 


3 


LattE 
CUBPACK 


0.0002 

0.0134 


0.0005 

0.0145 


0.0009 

0.0018 


0.0016 

0.0054 


0.0043 

0.0234 


0.0073 

0.0219 


0.0144 

0.0445 


0.0266 

0.0699 


0.0453 

0.1170 


0.0748 

0.2420 


4 


LattE 
CUBPACK 


0.0003 

0.0042 


0.0012 

0.0134 


0.0018 

0.0028 


0.0044 
0.0019 


0.0121 
0.0076 


0.0274 

0.5788 


0.0569 

4.7837 


0.1094 

4.3778 


0.2247 

22.3530 


0.4171 

54.3878 


5 


LattE 
CUBPACK 


0.0005 

0.0013 


0.0008 

0.0145 


0.0048 
0.0048 


0.0108 

0.0217 


0.0305 
0.0027 


0.0780 

37.0252 


0.0800 

128.2242 









Table 9: CUBPACK scaling with increased relative accuracy. "Relative Error" is a user- 
specified parameter of CUBPACK; "Expected Error" is an estimate of the absolute error, 
produced by CUBPACK's error estimators. Finally, the "Actual Error" is the difference of 
CUBPACK's result to the exact integral computed with LattE integrale. 



Relative Error Result Expected Error Actual Error # Evaluations Time (s) 



10" 


-2 


1260422511.762 


9185366.414 


94536.015 


4467 


0.00 


10" 


-3 


1260507955.807 


1173478.333 


9091.974 


9820 


0.01 


10" 


-4 


1260516650.281 


123541.490 


397.496 


34411 


0.04 


10" 


-5 


1260517042.311 


12588.455 


5.466 


104330 


0.10 


10" 


-6 


1260517047.653 


1257.553 


0.124 


357917 


0.31 


io- 


-7 


1260517047.691 


126.042 


0.086 


1344826 


1.16 


10" 


-8 


1260517047.775 


12.601 


0.002 


4707078 


4.15 


10" 


-9 


1260517047.777 


1.260 


< 10~ 3 


16224509 


14.09 


10" 


10 


1260517047.777 


0.126 


< 10~ 3 


55598639 


48.73 



which had every test be within 0.01 (the minimum time discrepancy recognized by its 
timer) of every other test case. 

CUBPACK differs from LattE integrale in that since it is based on numerical approxi- 
mations, one can ask for different levels of precision. Table 9 illustrates how CUBPACK scales 
with requested precision on a single, 4-dimensional, 10 linear form test case. It seems that 
CUBPACK scales linearly with the inverse of the requested precision — 10 times the precision 
requires about 3 times the work. All reported tests were done by expanding the multipli- 
cation of linear forms, and coding a Fortran 90 function to read in the resulting polynomial 
and evaluate it for specific points. 

5. One application: Voting theory 

Computation of integrals of polynomials over polyhedral regions is fundamental for 
many applications, including combinatorics, probability and statistics. In this last section 
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we wish to demonstrate the power of LattE integrale by attacking problems arising in the 
social sciences. In the mathematical theory of voting it was observed that the probability of 
events that can lead to singular election outcomes can be modeled as the number of lattice 
points inside a polytope divided by the number of lattice points of a simplex (see [46] and 
the references therein). Note that both the polytope and the simplex dilate proportional 
to the number n of voters. It is very well-known from the theory of Ehrhart functions that 
the counting functions are quasipolynomials (polynomials with periodic coefficients) that 
depend on n [22]. Thus when the quotient is evaluated the answer is asymptotically equal 
to the quotient of the leading coefficients of the two Ehrhart quasipolynomials involved. 

To illustrate this, consider the following example from [46]: There are three candidates 
a, b and c, and let the preference orders of the n = Yut=i Ui v °ters be 

abc (ni), acb (n 2 ), bac (n^), bca (714), cab (715), cba (no). 

Here, there are n\ voters who rank candidate a as first, b second, and c third, n 2 voters 
who rank b first, a second, c third, etc. Under simple plurality voting, the candidate with 
the most votes wins. But in a plurality runoff system, if no candidate wins more than 50% 
of the vote, the two candidates with the highest vote count advance to a second voting 
round. In [46] , the authors compute the probability that the simple plurality and plurality 
runoff systems give different winners. This requires setting up a system of equations that 
describes the situation that a wins by plurality but, using plurality runoff b obtains higher 
score than c and a majority of voters then prefer b to a. 

< ni + n 2 — n-s — n 4 

< n 3 + n 4 - n 5 - n 6 , 
1 

-- < -ni - n 2 - n 5 , 

1 = ni + n 2 + n 3 + n A + n 5 + n G , 
< Ui, i = 1, . . . , 6. 

This is done by computing the Ehrhart quasi-polynomial of the above polyhedron and 
dividing by the Ehrhart quasipolynomial of the simplex { (ni, n 2 , . . . , rig) : n\ + n 2 + • • • + 
n 6 = 1, rii > } (which is the space of all possible voting possibilities assuming that all 6 
rankings of three candidates are equally likely). All must be multiplied by 6 because the 
plurality winner may be a, b or c and the second position could be c not just b. As the 
authors observed, asymptotically, the leading coefficients of these two quasipolynomials is 
all that matter. In the concluding remarks the authors then posed the challenge of pushing 
the limit of such calculations for four-candidate elections which they observed is too big 
for their calculations. 

However, we have observed their calculation can be further simplified and accelerated 
because it is very well-known (see [22]) that the leading coefficient of the quasipolynomial 
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is always equal to the volume of the polytope with n = 1, thus one can directly perform 
the calculation of the volume (the volume of the simplex is well-known) and do a quotient 
of two numbers. The key step in finding the probabilities requires only finding the volume 
directly. Our algorithm corroborates that for the previous example the volume is 41 7 4720 , 
and when multiplied by 6 x 120 gives the probability these two voting systems give different 
winners for a large population: 12.33%. 

Using our code for exact integration we tackled the same problem for four candidates. 
In this case we have 24 variables associated to the orderings 



abcdini), abdc{ri2) , acbd{n^) , acdbin^), adbc(n 5 ), adcbijio) 
bacd(n 7 ), badc(n 8 ), bcad(n 9 ), bcda(n w ), bdac(nu), bdca{n\2) 
cabd(ni 3 ), cadb^riu), cbad(ni^), cbdain^), cdab(rin), cdba(rii 8 ) 



dabc(nig), dacb(ri2o), dbac(ri2i), dbca(ri22), dcab(n2s), dcba(n24) 

The equations and inequalities associated to the problem codify the following facts: 
The sum of all variables rii must be equal to the total number of voters. We have four 
inequalities expressing that when a is the plurality winner, b obtained a score higher than c, 
and c obtained a score higher than d, thus 

n\ + n 2 + n 3 + n A + n 5 + n e > n 7 + n 8 + n 9 + n w + n u + n 12 , 

n 7 + n 8 + n 9 + n w + n u + n 12 > n 13 + nu + ™i5 + n ie + n xl + n 18 , 

n 13 + n u + n 15 + n w + n 17 + n 18 > n 19 + n 20 + n 2 i + n 22 + n 23 + n 24 - 

These inequalities assume that the order was a > b > c > d but the answer we get 
should be multiplied by 4! = 24 to take into account other possible orders. Finally, we 
have to express the fact that but a majority of voters prefer b over a and that a did not 
achieve more than 50 percent of the vote (ni + n 2 + + n A + + < n/2). The volume 
of this polytope when n = 1 is 

2988379676768359 
7552997065814637134660504411827077120000' 
The probability is then the volume times 4! divided by the volume of the simplex 

{ (m, n 2 , • • • , n 24 ) : J2i n i = !> n i > }> 
which equals ^y. After a minute of computation using the cone decomposition method, 
we obtain the probability is 12.27%. 

We can continue the example by considering the same problem for five candidates. 
The five-candidate polytope has 5! = 120 variables. However, after LRS [25] enumerated 
over 12.5 million vertices, we terminated the program and decided the polytope is beyond 
our limits. We close by mentioning that after the first version of this paper was made 
available other authors proposed new ideas to compute these values using symmetries of 
the problem. See [47] . 
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